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Profnco 


At  present  very  little  information  is  available  about 
tho  effects  of  finite  length  on  the  spaco-chargo  limiting 
current  in  a drift  tube.  This  report  provides  accurate 
numerical  values  for  this  limiting  current  for  drift  tubo3 
of  several  longth3  and  for  oloctron  beams  of  several 
geometries  and  energies.  Those  values  may  bo  used  to  check 
various  analytical  approximations  for  tho  limiting  current 
and  in  the  do sign  of  experimental  equipment. 

Thanks  are  due  to  Profossor  Philip  Nielson,  Thomas 
denoni,  and  William  Proctor  for  their  guidance  on  this 
project  and  to  Frofossor  John  Jonos  for  his  assistance  with 
tho  numerical  technique. 
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Lint  of  Symbols 


tho  electrostatic  potential 

* 

the  eloctrostatic  potential  in  units  of  e/mcc 
1.6  x 10  1 ^coul. 
total  beam  current 

total  beam  current  in  units  of  e/m c^ 
the  space-charge  limiting  current 
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Abstract 


Tho  equation  governing  the  steady  state  potential 
distribution  crcntod  by  a cold  unneutralizod  annular 
relativistic  oloctron  boan  propagating  axially  ui thin  a 
finite  length  drift  tube  with  grounded  anode  foil,  collector 
plato,  and  walls  and  immersod  in  an  axially  directed  infinite 
magnetic  guide  fiold  is  solved  numerically  using  the  mothod 
of  finite  elements.  Those  numerical  results  aro  used  to 
check  tho  accuracy  of  available  analytical  approximations 
to  tho  limiting  current. 
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I Introduction 


background 

Many  applicr.tions  ol'  intense  re  inti  via  tic  electron 
boons  ( x-ray  and  microwave  generation*  ther  onucleor 
fusion  experiments,  and  collective  ion  acceleration 
schemes,  among  others)  roquiro  that  the  beam  bo  injected 
into  an  evacuated  drift  tube.  If  an  external  magnetic 
guide  field  is  applied,  the  injected  electrons  may 
propagate  thru  the  drift  tube  ns  a stable  beam.  If 
the  external  magnetic  field  is  sufficiently  strong 
(as  defined  in  Section  II)  the  maximum  current  which 
can  bo  totally  transmitted  thru  the  drift  tube  in 
steady  state  is  linitod  by  the  electron  number  density 
(space-charge)  at  the  center  of  the  tube.  For  currents 
above  the  s ace-chargc  limiting  current,  the  spac.o- 
chnrge  at  the  center  of  tho  tube  becomes  large  enough 
to  reflect  part  of  the  electron  beam.  In  this  papor 
the  effects  of  the  finite  length  of  tho  drift  tube 
on  the  space-charge  limiting  current  will  bo  considered 
numerically. 

For  a drift  tube  (Fig.  1)  with  L^Ii  and  with 
Bq  offectivly  infinite,  Bogdnnkovich  and  Rukhadzc 
(Ref.  1)  formulated  an  analytical  noproxination  for 
tho  spaco-charge  limiting  current  of  a solid  electron 
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Figure  1.  An  rlcctron  drift  tubo  of  finite  length  (L) 
in  an  effectively  infinite  axially  directed  magnetic 
guide  fi^ld  (B  )•  Tho  electrons  are  injected  nt  s=0 
and  occupy  tho°drift  tube  between  r=a  and  r“b„  (Ref.  2) 


bean  (Fig.  2a).  Their  a: nroxination  uns  derived  by 
interpolating  between  the  exact  solutions  for  nonrol ativistic 
and  ultrarelativi stic  olectron  beams.  Ilov/evor,  both  exact 
solutions  were  dorived  assuming  that  the  nunbor  density 
di 3trihution  shown  in  Fig.  2a  propagated  unchanged  thru 
the  drift  tube.  The3o  solutions  uero  derived  subject  to 
tho  further  assumption  that  L:»R  and,  thorefore,  the 
potential  is  constant  in  the  z direction  near  tho  tubo 
center  as  shown  in  Fig.  2c.  Proctor  and  Gononi  (Ref.  3) 
later  dorived  an  analytical  approximation  to  the  limiting 
current  analogous  to  tho  one  derived  by  Bogdnnkevich  and 
Rukhadse.  However,  in  their  derivation  annular  olectron 
beams  (Fig  2b)  woro  considered  in  addition  to  solid  beams. 


Figure  2.  Electron  number  density  distributions, 
o)  rndini  number  density  distribution  of  u solid  electron 
bean  at  the  injection  anode. 

b)  rndini  number  density  distribution  of  an  annular  (hollow) 
electron  beam  at  tho  injection  anode. 

c.)  axial  number  density  distribution  of  an  electron  beam 
propagating,  in  a long  (L  5R)  drift  tube  (Rof.  2). 

Furthermore,  their  derivation  does  not  require  that  the 
electron  number  density  propngnto  thru  the  drift  tube 
unchanged.  Tho  expression  for  the  limiting  current  in  n 
long  drift  tube  (as  derived  by  Proctor  and  Gononi ) is  given 
in  Appendix  A and  i:as  used  to  check  the  numerical  solutions 
of  tho  long  length  limit  of  tho  equation  derived  in  Section  II, 
Voronin  c_t  si . (Ref.  ) relaxed  the  restriction  that 
the  drift  tube  be  long  compared  to  tho  tubo  radius  and 
dorived  a rigorous  upper  bound  to  tho  limiting  currant 
of  a solid  cloctron  beam  filling  the  drift  space  (Fig  2a 
with  b=R).  Proctor  and  Gononi  (Rof.  5)  oxtondod  this 
analysis  to  include  a length  dependent  uppor  bound  to  tho 
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limiting  curront  of  annular*  electron  beams  (Fig.  2b), 

In  both  of  those  analyses  the  current  density  (nr)  and 
not  tho  olectron  number  density  (n)  was  assumed  constant 
in  z.  It  should  bo  emphasized  that  tho  Proctor  avid  Gononi 
formula  given  in  Appendix  A is  a rigorous  upper  bound 
but  not  necessarily  a least  upper  bound  to  tho  space-charge 
limiting  current.  It  is  not  an  expression  for  the  limiting 
current, 

Proctor  and  Gcnoni  (Ref.  6)  dcrivod  an  analytical 
approximation  for  tho  limiting  current  of  an  electron 
beam  with  an  electron  number  density  proportional  to  a 
delta  function  in  the  radial  direction  (Fig,  2b  with  a-»b). 

This  is  still,  howevor,  ossontinlly  a one  dimensional 
problem,  Millor  (Ref,  7)  has  developed  a fully  two  dimensional 
approximation  for  tho  limiting  curront  of  solid  electron 
beams  propagating  in  finite  length  drift  tubes  (given  in 
Appendix  A).  In  this  derivation,  however.  Miller  assumed 
that  tho  electron  number  density  remained  constant  throughout ■ 
the  drift  tubo,  Furthor,  ho  defined  the  limiting  curront 
as  that  current  which  produood  a potential  drop  (due  to 
the  space-charge)  botweon  tho  injection  anode  (z=0)  and 
the  tube  contor  (z=L/2)  equal  in  magnitude  to  the  olectron 
injection  energy.  Hovover,  it  has  boon  3hotm  (Ref,  2)  that 
tho  potential  drop  at  the  limiting  curront  nood  not  bo 
that  groat. 
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Objectives 


The  objectives  in  this  paper  are  to  present  numerical 
values  for  the  space-charge  limiting  currents  of  electron 
beams  propagating  in  finite  length  drift  tubos,  The 
limiting  currents  of  two  different  geometry  solid  beams 
with  eaoh  geometry  having  three  different  energies  will  be 
determined  numerically.  These  values  will  be  used  to 
estimate  the  accuracy  of  the  analytical  approximation  given 
in  Ref*  7*  Also,  two  different  hollow  beam  geometries 
with  each  geometry  having  three  different  energies  will 
be  considered.  The  limiting  currents  determined  for  these 
oases  will  be  presented  mostly  as  an  aid  to  future  development 

Approach 

la  Section  II  the  differential  equation  with  boundary 
oonditions  describing  the  electric  potential  inside  a 
finite  length  drift  tube  is  derived  subject  to  several 
simplifying  assumptions: 

1 ) The  drift  tube  and  electron  source  are  immersed 
in  a.  large  (as  defined  in  Section  II),  axially 
direoted  magnetic  guide  field* 

p 

2)  The  electron  kinetic  energy  ((/''-  1 )mc  ) is 
large  enough  (as  described  in  Section  II)  to 
make  scattering  from  the  anode  foil  at  z=0 
negligible* 

3)  The  electron  volocity  is  the  same  for  all 
electrons  at  z=0  and  is  totally  direoted  in 
the  +z  direction. 

4)  The  current  density  is  constant  aoross  the 
eleotron  beam  (r=a  to  r=b  in  Pig*  1 ) at 


Tho  numerical  analog  to  the  equation  derived  in  Sootion  IX 
is  dorivod  ir.  Section  III  ant)  tho  nurooricnl  method  (finite 
element  method)  used  to  solve  thin  equation  ia  discussed. 

Tho  limiting  currents  for  several  beam  geonotrios  and 
onergien  nro  given  in  Section  TV.  The no  limiting  eurronts 
are  comparod  to  tho  rigorous  upper  bound  for  the  limiting 
currents  given  in  Hof,  5»  l?or  the  solid  beams  tho  limiting 
currents  dotorm.lned  numerically  are  compared  to  tho  analytical 
approximation  given  in  Hof.  7.  Finally,  in  Soction  V 
recommendations  for  further  research  nro  made. 


II  Tho  Governing  Fount ion 


Dorlvrtl on 

Maxwell' 3 steady  state  equations  for  the  electro- 
magnetic fields  In  the  drift  tube  shown  in  Fig.  1 (in 
gaussian  units)  aro 


VXE  = 0 (In) 

V-E  = l^ven  (1b) 

VXH  = ^ nev  (1c) 

VS  =0  (Id) 


whore  n nnd  v aro  tho  electron  number  density  and  volocity. 
In  steady  stnto 


E = -V^  (?.) 

Equations  (1)  and  (?)  can  be  combined  to  givo  Foisson's 
equation 

Sftp  = liwenU(r)  (3) 

whore  U(r)  is  a unit  stop  function  such  that 


U(r) 


/ 0 a > r,  b<r 
\l  n < r s»b 


(U) 


By  electrically  connecting  tho  injection  anode,  collector 
plate,  and  drift  tube  walls  and  by  holding  then  all  at  a 


constant  potentinl  (dofined  to  be  zoro),  tho  boundary 
conditions  on  Eq  (3)  are 

0(r,O)  = 0(r,L)  = 0 (5ft) 

0(R,z)  = 0,  0(0,  s ) finite  (5b) 

Tho  divorgonce  of  Eq  (1c)  yiolds  the  equation  for  tho 
steady  state  conservation  of  charge : 

\7(nv)  =0  (6) 

If  tho  electrons  aro  injoctod  with  axially  directed  velocity, 
arc  not  scattered  significantly  by  t?io  nnodc  foil,  and  the 
drift  tubo  is  inmersod  in  a sufficiently  strong  magnetic 
guide  field,  Eq  (6)  reduces  to 

fr(nv)  = 0 (?) 

which  menn3  that  the  charge  density  (nv)  is  independent 
of  axial  position.  Theroforc,  sinco  nv  is  radially  constant 
at  the  injection  anode  (z=0),  it  is  constant  throughout 
tho  drift  tubo: 


nv  - n v = constant  (8) 

o o 

by  multiply 'uv  and  dividing  tho  right  hand  ride  (HUS)  of 
Eq  (3)  by  tho  velocity  of  an  electron  nnd  applying  Eq  (8), 
Foisson’s  equation  can  bo  expressed  in  toms  of  tho  current 
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density  (env  = enQvo)  and  the  electron  velocity  as 


V2^  = l ivr*  “0“0  (9) 

The  beam  current  (I)  is  computed  from  the  current  density 

by 


which  reduces  to 


b 


2 -IQ 

a 0 


env) (rdrd©) 


(10) 


I = en  v (b2-a2)  (11) 

c o 

for  cylindrical  geometry  and  constant  beam  area.  Combining 
Eq  (11 ) and  the  velocity  written  in  terms  of  the  relativity 
factor  (v  = cVw^  where  T - 1/V1  ” (v/c)2  ) Eq  (9) 
can  be  expressed  as 


il 


b2-.2  c 


(12) 


The  final  step  in  the  derivation  is  to  express  the 
relativity  factor  in  terms  of  the  electric  potential. 

To  do  this  the  total  derivitive  of  the  potential  is  expressed 
in  terms  of  the  electric  field: 


+ I-V0 

b -y.E 
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(13) 


Fy  talcing  tho  total  time  derivative  of  tho  total  relativistic 
onergy  of  an  electron 


dt 


= "^(c^P*  + rn^c'1-) 


.2  > 


(1!;) 


or 


£ = c2p  dP 

* dt  “ c i-  dt 


(1?) 


By  expressing  £ as  yhc‘  and  P ns  /mv,  Eq  ( 1 53,' ) reduces  to 


, 2 d/  dE 
me  rr  = v*-rr 
at  — dt 


(16) 


By  forming  tho  dot  product  of  v with  tho  Lcrents  forco 


dP 


s-St  = -<*•<£  + 

Eq  (16)  becomes  (noting  that  v-(vxB)  ~ 0) 


(17) 


<ur 

dt 


v • E 


(18) 


111c 


P>quationa  (13)  and  (18)  con  now  bo  combined  to  form  a 
constant  of  the  motion: 


dt 


(r-“V0  ) = 0 

nc‘ 


(19) 


By  expanding  the  total  derivative  in  Eq.  (19)  into  tho  sum 
of  partial  derivatives,  and  by  noting  that  this  i3  a steady 


1 0 


1 


w 


stnto  problem. 


v.V(  X-  -^0)  = 0 (20) 

me''  1 

Howovor,  since  vis  totally  diroctod  in  the  + z direction, 
Eq  (20)  reduces  to 


|r  O'-  -M)  = 0 (21) 

^ 1 mc'~ 

Therefore,  since  the  potential  at  z~0  is  zero,  and  sinco 
tho  term  in  parenthesis  in  Eq  (21 ) is  constant  in  z. 


(22) 


where 


Y0  is  the 

Finally,  using  Eq  (22 


electron  normalized  energy  a 
),  Eq  (12)  c.-in  be  expressed 


t 

as 


s=0 


h i)  U(r) 

b -n2  >fl-C  K +<^)"" 


(23) 


where 

l)  = ol/mc^ 

= e 0/mc. 

It  is  numerically  moro  convenient  to  oxpress  tho  boundary 
conditions  from  Eq  (5)  in  terms  of  tho  potential  at  z~0 
and  tho  derivative  of  the  potential  at  s-L/2  a3 
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(24a> 

(2lj.b) 


</>(r,0)  = 0,  f^(r,L/2)  = 0 

(/>  (R,z)  = 0,  (p(R,0)  finite 
which  arc  valid  for  a symetric  potential.  This  change 
in  the  boundary  conditions  permit  grontcr  accuracy  in  the 
numerical  solution  without  increasing  tho  number  of  mesh 
points  used. 

D1 scussion 

It  ha3  been  shown  that  not  all  values  of // permit 
physically  meaningful  solutions  to  Eq  (23)  (Hof.  ^ ) . 

That  is,  solutions  for  which  the  potential  is  nowhero 
positive  and  zero  only  at  tho  boundary  of  tho  drift  tube. 

The  assumptions  used  to  derive  Eq  (23)  were 

1)  Tho  externally  imposed,  axial  magnetic  guide 
field  (B  ) may  be  considered  .infinite, 

2)  Tho  kinetic  energy  of  the  electrons  at  the 
injection  anode  .is  large  enough  to  make 
scattering  from  tho  anode  foil  negligible c 

3)  The  electron  velocity  is  tho  same  for  all 
electrons  at  z=0  and  is  totally  directed  in 
tho  +z  direction. 

l|.)  Tho  current  density  is  constant  across  the 
electron  beam  (r=a.to  r~b  in  Fig,  1 ) at  z=0. 

From  a practical  standpoint,  the  first  two  of  these  assumptions 

are  the  most  difficult  to  obtain.  Thode  et_  al . (Ref.  0) 

has  quantifiod  tho  terms  "considered  infinite"  and  "negligible 

scattering,"  Their  findings  are  that,  for  a typical  injection 

anode  scattering  can  bo  neglected  if  the  total  particle 

1 2 


2 

energy  is  groator  than  2mc  . This  is  equivalent  to  an 

accelerating  field  of  0.0  HV  for  electrons.  AI30,  they 

show  that  the  magnetic  guido  field  may  be  considered  infinite 

if  the  ratio  of  the  cyclotron  frequence  to  tho  plasma 

frequency  (w  /w , a measure  of  tho  gyroradius  of  the  particle) 
c p 

at  .the  injection  anode  is  creator  than  0.  In  this  ratio 

the  cyclotron  frequency  is  defined  by  w = oB  /me  and  the 

c o 

plasma  frequoncy  is  defined  by  wp=  V . For  a typical 
electron  number  density  (10l2/cra-*)  this  requires  a magnetic 
guido  field  greater  than  about  70  kilo-gauss. 
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Ill  Finite  Element  Solution 


Derivation 

The  Unite  element  method  was  used  to  solve  Eq  (23). 

This  numerical  technique  involves  approximating  the  unknown 
solution  ( ) to  a differential  equation  with  an  interpolating 

function  (^?  ****): 


NM 


-l  l 


if*  (25) 

where  the  c^j’s  are  undetermined  constants  and  the  S^’s 
are  known  functions  (called  basis  functions)  satisfying 
the  same  boundary  conditions  as  (/  . The  method  chosen,  to 
determine  the  constants  was  Galerkin's  method  of  weighted 
residuals*  By  this  method  the  approximate  solution  to  the 
differential  equation  is  substituted  for  the  exact  solution, 
the  equation  is  then  multiplied  by  each  of  the  basis  functions 
in  turn,  and  each  equation  is  integrated  over  the  range  of 
variables*  One  such  equation  is  shown  below: 


R L/2 

/./!» 


3r^r 


+ d2vm^  Qij 


r*0  z=0 


■/ 


~d*c 

L/2 

U(r) 


) S1J  r dr  dz 


hi)  r dr  dz 


r«=0  z=0 


r1-(£  + ^il)“<: 


(26) 


This  produces  N times  M simultaneous  equations  for  the 
°ijts*  To  8°lye  these  equations  the  basis  functions  must 
be  specified.  For  the  present  problem  the  bilinear  Lagrange 
polynomial  basis  functions  given  by  Pr enter  (Ref,  9*128) 
were  used.  This  basis  is  defined  by 


pi*rr 

*“*.1-1 

rlS  r sri+1 

ri+i-ri 

Vlj-i 

*3-1S  X5*J 

ri4.l’P 

*1+1** 

pi  — p -ri+1 

ri+1-ri 

ZJS*-ZJ*1 

(27) 

P“ri-1 

pi-1^p^pi 

ri-*i-1 

'r'j-, 

zj-iSz^zj 

P~Pi-1 

pi-i-  r-ri 

ri-Pi-1 

ZJ-  z -z3+i 

0 

elsewhere 

One  element  of  this  basis  is  shown  in  Fig.  3«  These  linear 
polynomial  basis  functions  were  chosen  because  they  are  among 


Figure  3.  A canonical  bilinear  Lagrange  polynomial 
over  its  region  of  support  (Ref.9)« 
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the  simplest  to  ovalunto  numerically  and  because  they  produce 
a strongly  banded  system  of  algebraic  equations  to  be  solved. 
That  is  a 3yston  of  equations  in  which  most  of  the  coefficients 
aro  zero.  However,  using  this  basis  does  introduce  a difficulty 
with  the  second  derivative  - it  doesn't  exist.  This  difficulty 
can  bo  eliminated  by  integrating  the  loft  hand  side  (LHS) 
of  Eq  (26)  onco  by  parts: 


L/2 


lus  = 


z=0 


'( r u: 

1 d r 


in-i 


R 


R 


) dz  + 


r-0 


r=0 


(r  &T  ^ 
a z 


L/2 


) dr 


z=0 


n l/2 

r r , 2^j,m  i sw  Jsi. i 

J 1 T?  * 7?  fi  > ' 

r=0  z=0  ^ 

The  integrated  portions  of  Eq  (28)  can  bo  evaluated  using 
tho  boundary  conditions  for  the  differential  equation. 

■fiio  first  integral  is  obviously  zero  st  r equal  to  zero. 

For  r equal  to  R,  tho  basis  functions  must  all  be  zoro 
to  match  tho  boundary  conditions  o n </>.  Therefore,  tho  first, 
integral  is  identically  zero.  The  second  integral  also 
evaluates  to  zero  since  tho  basis  functions  arc  zero  at 
z oqual  to  zero  and  the  derivative  i3  zero  nt  z equal  to 
L/2.  Therefore,  tho  only  contribution  from  tho  LHS  of  Eq  (26) 
is  tho  third  integral  from  Eq  (28).  Combining  Eqa  (26) 


(28) 
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and  (28), 


R L/2 


- / / 


( 


r=0  2=0 


^_siJ 

a r a r 


a^NM  a 
^ z a z 


) r dr  dz 


R L/2 


U(r) 


l>  {)  v dr  dz 

bL'-a^  Vi-'u  +^i:iIr7" 


(29) 


f=0  5=0 

is  tho  actual  equation  to  bo  solved  numerically.  Since  i 
and  j vary  over  thoir  respective  ranges  of  values,  Eq  (29) 
actually  represents  a system  of  simultaneous  non-linear 
algebraic  equations.  This  set  of  equations  was  solved 
by  a non-linear  equation  solving  program  (Q1I)  written  by 
programmers  at  Sandia  Laboratories.  Each  of  the  equations 
represented  by  Eq  (29)  was  integrated  numerically  using 
assumed  values  for  the  c.  .'s  and  a two  dimensional  trapezoidal 
integration  approximation.  Those  equations  for  the  1 s 
were  solved  by  QN  using  a quasi-Newton  iterative  technique 
in  which  tho  equations  were  ropeatodly  linearized  and  solved 
until  successive  iterates  converged.  The  space-charge 
limiting  current  was  defined  numerically  as  the  largest 
current  for  which  the  iterative  equation  salvor  could 
produce  a converging  solution  that  was  physically  moangingf ul : 
that  is,  for  which  the  potential  waa  loss  than  zero  everywhere 
in  tho  interior  of  tho  drift  tube  and  for  which  the  potential 
root  the  boundary  conditions  specified. 
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Accuracy 

The  numerical  equation  (Eq  (29))  was  solved  usin^  a 
rao3h  siso  of  Il|  by  7 in  tho  r and  z directions.  The  mesh 
points  were  equally  spaced  in  the  z direction  and  equally 
spaced  in  tho  r direction  with  tho  exception  that  the 
nesh  points  falling  closest  to  the  inner  and  outor  boam 
radius  was  forced  tc  fall  exactly  at  that  radius.  To  check 
to  see  if  the  numerical  solution  was  converging  to  ono  value, 
tho  mesh  size  was  increased  up  to  the  limit  of  the  computer 
memory  (l£  by  15).  These  results  were  in  very  close  agreement 
uith  the  results  using  a llq  by  7 mesh.  The  limiting  currents 
were  computed  to  a precision  of  ^0.01  unite  of  normalized 
current  [)J  ) or  about  170  amps  of  electron  current.  In 
Table  I the  limiting  currents  found  by  solving  this  equation 
for  the  numerical  approximation  to  a delta  function  electron 
distribution  are  compared  to  the  results  given  in  Ref.  6. 

The  maximum  error  is  loss  than  2 percent.  Table  II  shov:s 
tho  limiting  currents  for  a drift  tube  of  length  £R.  At  this 
length,  the  limiting  curront  has  essentially  reached  its 
infinito  length  value.  Those  results  are  comoared  to  the 
infinite  length  results  given  in  Ref.  2.  Tho  maximum  error 
in  this  case  is  less  than  3 percent.  Based  on  those  comparisons 
the  maximum  error  in  the  numerical  results  for  finite  length 
drift  tubes  is  estimated  to  bo  on  the  order  of  5 percent. 
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Table  I.  Tho  limiting  curronts  for  a dolta  function 
distribution  of  olcotron3  is  compared  to  the  limiting  currents 
found  numerically  for  nn  electron  boam  having  a geometry  of 
a=0*4999,  b=0.5001  , Both  bonms  nro  qontered  about  r=0.5 
end  have  an  injection  onorgy  of  20mc. 


Percent 


Length 

Delta  Emotion 

Numerical 

Diff  oi*oi 

1 

7.30 

7.43 

1 .7 

2 

4 .60 

4. 62 

0.4 

3 

4.07 

4 .06 

0.2 

5 

3.01 

3.01 

0.0 

oo 

3.78 

3.75 

0.8 
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Tablo  II.  Tho  limiting  currents  fox*  electron  beams  of 
several  energies  and  gcometrios  propagating  in  nn  infinito 
length  drift  tube  are  compared  to  the  limiting  currents  for 
electron  boam3  propagating  in  a drift  tube  of  length  £R. 

A length  of  £R  was  used  since  that  was  the  longest  length 
in  which  any  offocts  of  finite  length  wox-o  obsorvod. 


Energy 

Limiting 

Current 

Percent 

(infinito 

) 5R 

Difference 

n=0.0,  b~0.5 

2.0 

0.22 

0.22 

0.0 

5.0 

1 .27 

1.29 

1.5 

20.0 

7.31 

7.37 

0.8 

n=0.0,  b=1.0 

2.0 

0.55 

0.56 

1.8 

5.0 

3.12 

3.16 

1.3 

20.0 

17.50 

17.75 

1 .0 

a^O.57,  b=0.60 

2.0 

0.43 

0.1(4 

2.1 

5.0 

2.53 

2.56 

1.1 

20.0 

15.19 

15.50 

2.0 

a=0. 30,  b-0.60 

2.0 

0.32 

0.33 

3.0 

5.0 

1 .67 

1 .91 

2.0 

20.0 

10.91 

11.12 

1.9 

IV  Results 


Annular  Beans 

Two  annular  boon  goomotries  wcro  considered.  They 
wore  tho  numerical  approximation  to  a delta  function 
distribution  of  electrons  (n^O.4999,  b=0.5001 ) and  a thicker 
beam  (a=0.4>  b=0.6).  Both  bonms  were  centered  about  r=0.5» 
Tablos  III  (delta  function)  and  IV  (thicker  beam)  givo  tho 
results  of  solving  Eq  (?9)  for  tho  opace-chargo  limiting 
current.  Also  prosentod  nro  tho  upper  bounds  given  by  tho 
formula  derivod  in  Ref.  5 and  reproduced  in  Appendix  A. 

It  should  bo  notod  that  in  all  cases  tho  limiting  current 
does  approach  the  numerically  correct  infinite  drift  tube 
limit  as  the  tube  length  increases  and  that  tho  rosults 
for  tho  thin  annulus  do  agree  with  provious  results. 

At  present  thorn  i3  no  analytical  approximating  fornmla 
available  for  tho  limiting  curront  in  an  annular  beam. 
Therefore,  tho  data  in  these  tables  nro  prosentod  S3  an  aid 
to  futuro  researchers. 


Table  III.  Limiting  currents  for  a geometry  of  a=0.4999» 
b=0.5001  (a  numerical  approximation  to  a dolta  function 
distribution  of  electrons). 


Length  Upper  Bound  Numerical 


1 

X =2.0 

0.71 

0.66 

2 

o.i|i+ 

0.41 

3 

0.38 

0.35 

5 

0.34 

0.33 

CO 

0.33 

0.32 

1 

X=8.0 

8.13 

7.43 

2 

5. 02 

4.62 

3 

4.33 

4.06 

5 

3.96 

3.81 

Ad 

3.75 

3.75 

1 

^ =20.0 

25.1 8 

22.61 

2 

H5.52 

14.10 

3 

13.39 

12.43 

5 

12.25 

11.76 

c<D 

11.59 

11.59 

CM  r'YlA  g t-  CM  fVlA  g CM  nMA  ^ 


Tablo  IV.  Limiting  currents  for  thicic  annular  beam3  (a=0.4, 
b-0.6 ) 


Length 


Upper  Bound 


X =2.0 


X =8.0 


X =20.0 
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Numerical 


0.86 

0.80 

0.50 

0.46 

0.42 

0.39 

0.38 

O.36 

0.36 

0.35 

9.93 

8.79 

5.72 

5.09 

I)..  06 

4.41 

4.41 

4.13 

4. 16 

4.02 

30.70 

26.43 

17.69 

15.32 

15.04 

13.31 

13.65 

12.53 

12.85 

12.28 

Solid  Boans 


Two  solid  beam  geometries  wore  also  considered  (filling 
and  half  filling  the  drift  tube)  and  the  results  of  solving 
Eq  (29)  for  these  geomotrios  are  given  in  Tables  V (filling) 
and  VI  (half  filling).  A3  with  annular  beams,  the  upper 
bounds  to  the  limiting  currents  are  also  presented  in  these 
tables.  However,  Hiller  has  derived  an  analytical  approximation 
to  the  limiting  currents  of  solid  beams: 


4<L>  - b 


where  3s  i3  the  n'th  zero  of  tbo  J Bo  seel  function, 
n o 

This  equation  is  approximated  by  Miller  for  tubes  having 
a length  greater  than  the  radius  by 


-1 


(30) 


^l(L)  = ^l(1-2o“  'Vj/2)  -1  (31) 

where  ^ is  the  infinite  length  limiting  current. 

Figures  4-9  show  Eq  (30)  as  compared  to  the  numerically 
dorived  limiting  currents  and  Figures  show  Eq  (31  ) 

as  compared  to  the  numerically  dorived  limiting  currents. 

It  should  bo  noted  that  tho  inf ini to  length  limiting  current 
in  Eq  (31  ) was  obtained  from  tho  Proctor  and  tienoni  formula 
(Ref.  3)  and  not  from  Minor’s  expression  for  that  current. 
Comparison  of  those  two  sots  of  graphs  3howa  that  Eq  (31  ) 
is  in  much  bettor  agreemoiit  with  tho  numerical  rosults. 
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Table  V. 

Limiting  currents  for 

n solid  beam  filling  the 

drift  tube 

. (a=0.0,  b=1 .0 ) . 

Length 

Upper  Bound 

K=2.0 

Numerical 

1 

1.76 

1!  .42 

2 

0.93 

0.69 

3 

0.75 

0.60 

5 

0.70 

0.57 

QO 

0.65 

£=5.0 

0.57 

1 

10. 44 

7.72 

2 

5.51 

3.76 

3 

4.59 

3.26 

5 

4.12 

3.10 

oo 

3.86 

<£  = 20.0 

3.08 

1 

63.80 

42.90 

2 

33.15 

20.89 

3 

27.64 

18.10 

5 

24.82 

17.21 

aO 

23.23 

17.12 
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Table  VI.  Limiting  currants  for  n solid  bean  half  filling 
tho  drift  tube.  (a-0.0,  b=0.5>) 


Numerical 


0.53 

0.26 

0.23 

0.21 

0.21 


3.10 
1 .51 
1 .31 
1.21; 
1.21; 


17.80 

8.67 

7.51 

7.14 

7.10 


V Conclusion 


Tho  results  shown  in  Section  IV  and  the  program  that 
generated  then  form  a comploto  numerical  colution  to  tho 
problem  of  determining  tho  space-charge  limiting  current 
in  a drift  tube  of  finite  length  imnorsed  in  a strong 
guide  field.  Those  rosults  show  good  agreement  with  the 
formula  dorived  by  Miller  in  Ref.  7.  However,  by  comparing 
the  results  from  tho  simpler  formula  (Eq  ( 31 ) ) using  the 
Proctor  and  Genoni  expression  for  tho  infinito  length 
limiting  current  (Rof.  3)  with  the  results  using  Millor’s 
expression  as  dorived  it  appears  obvious  that  tho  assumptions 
Miller  used  (n  constant,  at  tho  limiting  current  ]^/>  (0,L/2)f-  ^-1  ) 
and  Froctor  and  Genoni  did  not  use  are  less  accurate  then  the 
assumptions  used  in  thi3  papor  (nv  constant,  y?(P,L/2)|  — ¥a  -1  ) 

It  seoms  reasonable  that  any  extension  of  Miller’s  formula 
should  bo  based  on  tho  assumption  that  nv  and  not  n i3  constant 
in  z and  should  uso  a difforent  definition  of  the  conditions 
on  tho  potential  at  tho  limiting  current. 
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Appendix  A:  Analytical  Formula  from  tbo  Literature 

Tho  Upper  Bound  Formula  (Hof,  s) 

The  upper  bound  to  the  limiting  current  ( t)^)  is 
g.ivon  by 

t\TB(L)  0'2+  ^2  )(r//3-1  )3/2  (A1  ) 

with  k determined  by  the  boundary  conditions 

a2[jQ(ka)  + CYo(lca)]  + a^af^Oca)  -t  CY^ka)]  = 0 (A2a) 

b2  [J0(kb)  + CYo0;b)]  + b1  kb  [j^  ( kb ) + CYj(kb)]  = 0 (A2b) 

where  J and  Y are  Bossol  functions  of  the  first  and  second 
m m 

kind  of  ordor  n,  a and  b aro  tho  innor  and  outer  beam 
radius,  and  C is  a constant  to  be  determined  along  with  k. 
Tho  a^  and  b^  aro  given  by 


n1  = I0('rra/L^  (A3a) 

a2  = (fra/L)^  (IT a/L)  (A3b) 

b1  = KQ(ir/L)lo(frb  /L)  - Iq  ( 7r/L ) Kq  ( TTh/L ) (A3c) 

b2  = ( dr b/L ) Jkq {ir/L ) ( -Dr b/L ) + I^rr/L)^  (-TTb/L)]  (A3d) 


v;hore  I and  K ai'e  modified  Bessel  functions  of  the  first 
m n 

and  second  kind  of  order  m. 


41 


Appendix  B:  Finite  l'l  emont  Commit  or  Program 


Tho  program  that  solves  Eq  (29)  was  written  ns  a 

series  of  subprograms.  Those  subprograms  are  called  by 

tho  supervising  program  ( DRIVER ).  DRIVER  is  responsible 

for  reading  all  parameter  cards,  insuring  that  the  parameters 

are  valid,  computing  the  locations  of  tho  mosh  points, 

calling  the  appropriate  oquation  solver  (ONED  or  TY.'OD), 

and  printing  the  results.  The  valid  parameters  are 

RO  inner  radius  of  the  electron  beam 

RN  outor  radius  of  tho  electron  beam 

NR  numbor  of  mesh  points  in  the  r direction 

V beam  current.  This  is  an  optional  parameter 

Which,  if  omitted,  causes  DRIVER  to  solve 
for  the  limiting  current 
GAI'O  electron  injection  energy 
Z 31  length  of  the  drift  tube  (required  for 

two  dinonsional  runs) 

NZ  number  of  mesh  roints  in  tho  z direction 

(also  required  for  t wo  dimensional  runs) 

ID  1 or  2 for  one  or  two  dimensional  runs 

Subroutine  TWOD  (ONED)  is  called  by  DRIVER  to  solve 
two  dimensional  problem  (Eq  (29))  (one  dimensional  problem 
CEq  (29)  with  ^ z = 0)).  Tho  problem  is  solved  using  QN, 
a non-linear  equation  solver  writton  bv  tho  programmers 
at  Sendia  Laboratories.  QN  solves  a series  of  non-linear 
algebraic  equations  cf  the  form 


Fi(c1,c;„...ci)  = 0 ( B1  ) 

for  tho  vnluos  of  tho  c^.  Tho  subroutines  FOFRZ  (FOFR) 


1*3 


celled  by  QN  to  evaluate  these  non-linear  equations  which 
it  does  by  directly  integrating  Eq  (29)  using  subroutines 
RZINTn  (RINTn).  If  QN  is  able  to  produce  a set  of  c^ 
that  make  all  F^  small  (less  than  ERRMAX  or  10”  ^ in  this 
case),  QN  returns  to  TWOD  (ONED)  with  IFLAG  = 2,  otherwise, 

IPLAG  is  sot  to  some  number  greator  than  2.  If,  on  returning 
to  TUOD  (ONED)  IFLAG  = 2 the  eloctron  current  ( ) Is  increased 

( was  loss  than  the  limiting  current),  otherwise,  the 
current  i3  decreased  by  H (in  both  cases),  the  stop  size. 

II  is  then  halved  and  this  cycle  continues  untill  H is  les3 
than  0.01  (0.001  ) when  control  is  returned  to  DRIVER  and 
the  next  parameter  cord  is  processed. 
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PROGRAM  PR  IV;:  R ( INPUT.  OUTPUT) 

COmON/DATA/'RO.  RN,  ZN,  NR,  NRL 1 , NRL2.NZ,  V,  G0M3 
V C0ITOM/C0EF/CC201) 

COmOM/RRDlUS/RaOl) 

COir.OU/Lb't!GTH/Z(  101) 

REAL  K 

NAMEL I ST/VAR IBL/RO,  RN, ZN,  NR. HZ, V, GAM3, ID 
R0=G.0 
RN--0.0 
ZN=R . 0 
Gfttt3=0.6 
I CONTINUE 
V=0.0 
1D-0 

READ  VAR  I PL 

IF  (EOF (5L INPUT) .NE.0)  GO  TO  509 
IF  ( ( ID. LT. 1 ) .OR. ( 1D.GT.2) ) GO  10  900 
IF  (RN.LE.R0)  GO  TO  910 
IF  (IP. Ed. 1)  GO  TO  10 

IF  ( (NR&NZ.EQ. C) .OR. (NR*NZ.GT. 109) ) GO  TO  920 
10  CONTINUE 

IF  (NR.GT.1G0)  GO  TO  330 

GET  DATE  AND  TIME 

A=0.0 
0=0.0 

CALL  DATE (A) 

COLL  TIME  CD) 

COMPUTE  THE  INDEX  AT  THE  INNER  BEAM  EDGE 

NRL  1 = IF  I XCC . 9999-:  R0*FLOAT(NR) ) 4 1 
IF  CNRL1.LT. 2)  NRL! =2 
IF  (R6.E0.G.0)  NRL 1=1 

COMPUTE  THE  INDEX  AT  THE  OUTER  BEAM  EDGE 

NRL2" IFIX(RN*FLOAT CNR) ) 41 
IF  (HRL2.LE.NRL 1)  NRL 2= NRL 1+1 

COMPUTE  THE  R VALUES 

H--1 .0/FLOATCNR) 

NREND--NR41 
DO  100  I = 1,1  IRENE 
RC 1)  -FLOAT ( 1-1 MH 
ICO  CONTINUE 

FORCE  THE  INNER  AMD  OUTER  EDGES  TO  BE  RIGHT 

R(l!RL!)=RO 
RCNRL2) =RN 
IFLAG=0 

CALL  ERXSET (0, 0) 

IF  (1D.E0. 1)  GO  TO  200 


kS 


( 


c 

c 


150 


C 

C 

c 

200 


C 

C 

c 


300 


c 

c 

c 

c 

400 

uioo 


c 

c 

c 

500 

2U00 


5000 


3000 

550 


C 

C 

C 


COO 

4000 


COMPUTE  THE  VALUES  OE  Z AT  MUCH  UE  WANT  GAMMA  VALUES 
NZEND-NZ+1 

K"ZN/(2.0+EL0AT(NZ)) 
no  150  J 1 , NZF.ND 
Z(J)**f  LOATCJ-IHKK 
CONTI DUE 
GO  TO  300 

F I HP  THE  ONE  DIMENSIONAL  LIMITING  CURRENT 

CONTINUE 

CALL  ONI  P< I FLAG) 

C(NR+1 ) rGAITQ 
GO  TO  400 

FIND  THE  TUQ  D1IF.NS10NAL  CURRENT 

CONTINUE 

CALL  TUMP (1 FLAG) 

GO  TO  400 

FRINT  THE  RESULTS 

CONTINUE 

PRINT  1000, 1FLAG 
I URMA1  (1II1.6H1I  LAG-,  12) 

PRINT  5000 

IF  C ID.EU.2)  GO  TO  600 


PR  I NT  NIL  ONT  P II  TENS  IONA!  RESULTS 


CONTINUE 

PRINT  2000, V,GAMA,RO,RN,  A,0 

FORMAT i III  , "llll“",l‘10.3,5X,  "GAIT3"",F?’.2,5X,  "R0«',F6.4,5X,  "RN»", 
1I-6.4.5X, "DATE-  ",ftI0,5X, "TIME* ",HIU) 

IF  ( 1FL0G.GT.3)  GO  TO  1 
PRINT  5000 
FORMAT  UNO) 

DO  550  J - 1 , NREND 

PRINT  3000, J,RCJ) ,C(J) ,C( J) -CANO 

FORITATv  III  "I-’’.  15, 5X,  "R*  ",FC.4,5X,  "GAM  ■ " , I'l 6 . ?, DX,  "PHI»",F16.7> 
CONI  INAL 
GO  TO  ) 


PRINT  HIE  THU  DIMENSIONAL  RESULTS 


CONTINUE 

PRINT  4000 , V, GAMO, RU.RN, ZN, A . P. 

FORMAT ( 111  , "Nil-", I 5.3,5X,  "GAM  I ",I:7.7,5X. 
1E6.4,5X,  "LI  HA  I II  * ,F5.4,5X,  "DAll>  ".A  10,  NX. 
IF  (IFLAn.GT. 3)  GO  TO  1 


"R0"",r6.4,5X, 
"TIME"",  A 10) 


"RN* 


PRINT  5000 


do  650  J"1,h;:end 

DO  640  I - 1 , Nk'i'NP 
CP “G AMO 


IF  ( ( 1 . Hi: . NRi-.ND) . AND.  ( J . NE . 1 ) ) CP-C(0  I+f  l-D'i-NZ) 


J|6 


c 


PRINT  5000, 1.  J. R (I ), Z (J), CP, CP-GAM3 

5U0P  FORMATC 1H  , "I““, I5,5X,  "J  = “,  15, 5X,  "R=\F6.4,5X,  ,Z=\F9.G,5X,  "GAM”", 
1F15.7.5X, "PHI*“,F15.7) 

640  CONTINUE 
PRINT  9200 
650  CONTINUE 
GO  TO  1 
900  CONTINUE 

DIMENSION  NOT  ViTLID 

PRINT  9001 

9001  FDRHATC INI, " ID  NOT  SPECIFIED  OR  NOT  VALID.  ID  MUST  BE  EITHER  1 OR 
12.  RUN  TERMINATED. ") 

GO  TO  999 


C 

C 

c 


c 

c 

c 


c 

c 

c 


ILLEGAL  RO,  RN  VALUES 

910  CONTINUE 

PRINT  9002, RO. PM 

9002  FORMATC 1111, "STARTING  RADIUS  MUST  GE  LESS  THAN  THE  ENDING  RADIUS’, 
15X, "R0=”,F7.4,5X, “RN- \F7.4) 

GO  TO  999 

ERROR  IN  THE  NUM3ER  OF  POINTS  REQUESTED 

920  CONTINUE 

PRINT  9003, NR, N2 

9003  FORMATC 1H1, "THE  PRODUCT  OF  NR  AND  NZ  MUST  PE  LESS  THAN  100'. 

1"  AND  GREATER  THAN  0.  IF  MORE  POINTS  ARE  REQUIRED,  CHANGES’, 

2"  MUST  HE  MADE  TO  TUQD",/, 

220X, "NR”", 13, 5X, “NZ=", 13) 

GO  TO  999 

TOO  MANY  R VALUES  REQUESTED 

930  CONTINUE 

PRINT  9004, NR 

9034  FORMATC 1 HI, "THE  NUMBER  OF  R VALUES  REQUESTED  MUST  DE  LF.SS  THAN  IPO 
1.  IF  MORE  VALUES  ARE  REQUIRED  CHANGED  MUST  PE  MADE  TO  GHED", 

2SX, "NR”", 14) 

999  CONTINUE 
STOP 
END 


C 

c 

c 
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SUBROUTINE  01-0 ( IFLAG) 

EXTERNAL  FOFR 

COMMON/DATiVRO,  PM . ZN, NR,  URL  1 . NRL2,  HZ,  V,  Gftl'10 
comoH/coEF/r(  ltrzo) 

COmON^RDDIUS/RUOn 

C0MM0N/0NU3RK/D  ISMAXC 100) , WORK  (2 1060) , IUORKC 120) 
CONSTANTS  FOR  QH 


RELERR" 1E-? 

ABSERR=0.0 
RES-0.0 
1-13  AND  =1 
VLA5T-0.0 
IFLG2c0 

IF  (V. EU. 0.0)  GO  TO  BOO 

JEND= 1 

IFLG2“ 1 

HV=0.0 

GO  TO  050 

FIND  THE  I1AX  CURRENT 

830  CONTINUE 
HV=GAM0 
DO  855  I - 1 , NR 
C( 1) =GAMT 

D ISMAXC I )“GAM0-1  .0 
855  CONTINUE 
8G0  CONTINUE 
VLAST=V 
V-'V+HV 
HV=HV"!<2 . 0 

CALL  QN (FOFR, NR, C, M3AND, D ISMAX, RELERR, RDSERR, IFLAG, RES, 
1L0RK, I WORK) 

IF  ( IFLAG. LC. 3)  GO  TO  GCO 


COMPUTE  THE  NUMBER  OF  TIMES  THE  BINARY  SEARCH  MUST  DE  EXICUTED 
TO  GET  C.U01  ACCURACY.  THEN  ADD  1. 

TEND- I F 1X( ALOG C (V-VLAST) /(2 . 0*0 .001)) /ALOG (2 . 0) ) +1 
050  CONTINUE 
HVeV 
V=2.0*V 
IFLAG“7 

BINARY  SEARCH  THE  CURRENT.  IE  PN  RETURNS  A DID-NOT-CONVERGE 
RETURN  CODE  ASSUME  IIIE  CURRENT  IS  100  NIGH.  ELSE  THE  CURRENT  IS 
TOU  LOU 

DO  780  Jcl, JFNP 
DO  300  1-1, NR 
D ISMAX(  I ) •'-GAMy- 1 . 0 
300  CONTINUE 

IF C I FLAG. LE. 3)  GO  TO  500 
DO  100  I ■*  1 , NR 
C( I) rGAM3 


k* 


o n n n n n n 


100  CONTINUE 
500  CONTINUE 


FOR  THE  FURPOSE  OF  DOING  THE  BINARY  SEARCH  A REDUCTION  OF 
10  ORDERS  OF  MAGNITUDE  IS  CLOSE  ENOUGH 


IF( IFLAG.LE.3)  V-V+HV 
IF(  IFLAG.GT.3)  V=V-HV 
HV-HV/2.0 
IFLAG= 1 
400  CONTINUE 

CALL  QN ( FOFR, HR , C , M3AND, D I SMHX. RELERR, ABSERR, IFLAG.RES, 
1U0RK,  ILJORIO 
IF ( I FLAG. LE. 31  VLAST=V 
200  CONTINUE 

FOR  THE  FINAL  PASS  THE  ERROR  MUST  DE  AS  SMALL  AS  POSSIBLE 

IF ( IFLAG.LT. 3)  GO  TO  900 
IF ( 1FLAG.E0.3)  GO  TO  700 
IF  C lFLGZ.CCJ . 1 ) CO  TO  900 
V-VLflST 
DO  GOO  1=1, HR 
CCD  =GAM3 

D ISMAXC I ) =GAM3” 1 . 0 
600  CONTINUE 
IFLAG-1 
700  CONTINUE 

CALL  ON (FOFR, NR, C , MO AND. D ISMAX, RELERR, ABSERR, IFLAG.RES, 
lldORK,  IUORIO 

IF  C IFLAG.E0.3)  GO  TO  700 
see  CONTINUE 
RETURN 
END 
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( SUBROUTINE  FOFR(C.F) 

COMMON/DATR/RO , RH , ZN , NR , NRL  l ; NRL2,  NZ , V,  GAM0 
COMMON/RAD 1US/R  (101) 

COMMON/SUB/ I FLAG, M, G (?) 

DIMENSION  F (5000). C (5000) 

M=7 

DO  000  1=1. NR 

INSURE  THAT  ALL  C(I)  ARE  GREATER  THAN  1.0 

IF(C( I) .LE. l.O)  GO  TO  900 

INTEGRATE  THE  SECOND  HALF  INTERVAL 

F( I) =RINT2( I.C) 

IF  ( 1FLAG.NE .0)  GO  TO  900 

IF  R IS  GREATER  THAN  0.0  INTEGRATE  THE  FIRST  HALF  OF  THE 
INTERVAL 

IF  (l.GT.l)  F( I)  =F( D+RINT1  ( I.C) 

IF  (IFLAG.NE.O)  GO  TO  900 
BOB  CONTINUE 
GO  TO  999 

ASSUME  THAT  V IS  GREATER  THAN  V(l I MIT) 

AND  FORCE  ON  TO  CALL  THE  SERIES  DIVERGING 

< 900  CONTINUE 

DO  901  1=1, NR 
F(I)=1 .0 
901  CONTINUE 
999  CONTINUE 
RETURN 
END 
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FUNCTION  RINTl(I.C) 

common/pat  a/ru  , rn  . zn, hr, nrl i , nrl2  , nz  , v,  cnra 
common/rad  i us/r  non 

COMMON/SUB/IFLAG,M,G(?) 

DIMENSION  CC-JOl) 

H"R(  I)-R(  1-1) 

THE  DER1VITIVES  ARE  CONSTANT  OVER  EACH  SUB  INTERVAL 

PRI“1 .0/(R(  I)-R(  I-D) 

DR11—  DR  I 

COMPUTE  THE  FUNCTION  VALUES  AT  EACH  POINT  IN  THE  SUB  INTERVAL 
FOR  THE  SIMPSONS  RULE  INTEGRAL 

DO  10U  J»1,M 

RA'R  ( I - 1 ) -.'-(FlOATC  J- I )*!()/(  FLOATdl- 1 ) J 
G ( J ) » CC  ( I - 1 ) *DR  I H C ( I ) *DR  I ) *DR I * KA 

IF  ME  ARE  IN  THE  PART  OF  THE  TUBE  OCCUF IFD  BY  THE  BEAM 
ADD  IN  THE  NON-LI HEAR  TERM 

IF  ( ( I .LE.NRL1) .OR. ( 1 .GT.NPL2) ) GO  TO  100 
Rl"(Rfi-R<I-l))/(R(I)-R(I-l)) 

R 1 1 = (R( J ) -RA) /(R (I)-R(I-!)) 

GAM-c(i-nxrv,ri+c(n>!Ri 
IF  (GAM. LE. 1.0)  GO  TO  BOO 

G(J)  -G(J)-'-(4.  B*V*GAM  :<R  1 *RA )/(( SORT  ( GAM**2- 1 . 0) ) * ( RN**2-R0*#2 ) ) 
100  CONTINUE 

DO  A S I UPSON  RULE  INTEGRAL 

R1HT1 =G( 1 ) 

MEND-M- 1 

DU  200  J “2, MEND, 2 
R INTI  *=R  JNT1+4.  U*G(  J) T-2.0>KG(J+1 ) 

200  CONTINUE 

RIHT1  '--R1NT1-GCI1) 

RINT1  nRiNTl*(H/(3.3*FL0AT(M~l) ) ) 

I FLAG =0 
GO  TO  399 

NEGATIVE  SQUARE  RUOT///NOTIFY  CALLING  PROGRAM 

900  CONTINUE 
R INTI -0.0 
I FLAG- 1 
999  CONTINUE 
RETURN 
END 
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FUNCTION  F;INT2(I.C) 

common/d ata/ro , rn . zn  , nr  , nrli  , nrl2  , mz  , v,  g ai  ;o 

COMMON/RAD IUS/R ( 101) 

COHMON/SUB/ IFLAG.  H,  G(7) 

D I TENSION  C (20 1 ) 

H°R( I+l)-R( I) 

THE  DER IVIT1VES  ORE  CONSTANT  OVER  EACH  SUB  INTERVAL 

DR1=- 1 .0/(R(  I + D-RC  I) ) 

DRI1— PR1 

COMPUTE  THE  FUNCTION  VALUES  AT  EACH  POINT  IN  THE  SUB  INTERVAL 
FOR  THE  SIMPSONS  RULE  INTEGRAL 

DO  100  J = 1 , M 

RA'-RC  I)  -K  FLOAT  (J--  1)*H)  ••'(FLOAT  (M-l) ) 

G ( J ) = (C  ( I ) Ji'DR  I +G AMQ*DR  1 1 ) -DR  I H;RA 

IF  (I.LT.NR)  G(J)  = (C(I):;T>RI+C(1H)*DRI1)*DRI*RA 

IF  UE  ARE  IN  THE  TART  OF  THE  TUBE  OCCUPIED  BY  THE  BEAM 
ADD  IN  THE  NON-LINEAR  TERM 

IF  ( ( I .LT.NRL1)  .OR.  (I  .GE.NRL.2) ) GO  TO  100 
R1  = (R(  H I) -RA) /(R(I  + 1)-R(D) 

R 1 1 " (Rfl-R  (1))/(R(1  + 1)-R(D) 

GAM-C  ( I ) *R  1 +GAM3TR 1 1 

IF  (I.LT.NR)  GAii:  0 ( I) ’KR I+C ( 1+1 ) *R  1 1 

IF  (GAM. LE. 1.0)  GO  TO  BOO 

G(  J)  -G(J)  H4.Q*V*GAI1*RI*RA)/(  CSGRT(GAN**2- 1 . 0) ) *(RN**2-R0**2) ) 
100  CONTINUE 

DO  A SIMPSON  RULE  INTEGRAL 


RINT2~G( 1) 

MEHDeM- 1 

DO  230  J -2, MEND. 2 
R I NT2 -R 1 HT2+4 . 0*G (J)+2.0*G(J+1) 

200  CONTINUE 

RINT2“R INT2-G (M) 

R INT2“R  INT2*(H/(3.0>i!FL0AT(M- 1 ) ) ) 

1FLAG=0 
GO  TO  BSD 

NEGATIVE  SQUARE  R00T///N3TIFY  CALLING  PROGRAM 

900  CONTINUE 
RINT2C0 . 0 
I FLAG -1 
9B9  CONTINUE 
RETURN 
END 
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SUBROUTINE  TUOD(IFLAG) 

EXTERNAL  FOFRZ 

C OMMQ  N y'D  ft  Ti VR  0 , RN,  ZN,  NR,  NRL 1 , NRL2,  NZ,  V,  GAM3 
COMMOiVRADIUS/R(  1 1) 

COmON/LEHGTH/Z  (11) 

COMMON/COEF/COOl) 

COmON/QHUORK/D ismxc  100)  ,U0R;<(210e0) , ItJ0RK(  12(0 
IFLG2k0 

IF  (V. EQ .0.0)  IFLG2= 1 
PRINT  9399 
9999  FORMAT C 1H1) 

VLA5T-V 

HV=GAM0 

CONSTANTS  FOR  QU 

ITBAND5!MZ+1 
RELERR=lE-7 
ABSERR=U. 0 
RES-0 . 0 

CHECK  FOR  ONE  PASS  PNLY 
IFLAG=7 

IF  ( IFLG2.EO.0)  GO  TO  GOO 

COMPUTE  THE  INITIAL  TWO  DIMENSIONAL  GUESS 

50  CONTINUE 
NEND=NR*NZ 
DO  100  IK1,NEND 
C( I) =GAM0 
100  CONTINUE 

FIND  1HE  MAX  CURRENT 

150  CONTINUE 

DO  203  NM.NEND 
DISMAX(N) =GAM3-1 .0 
200  CONTINUE 
IFLAG“3 
VLAST=V 
V'=V+HV 
HV=HV*2 . 0 

CALL  ON (FOFRZ, NR:«Z,C,MBAND, DISMAX,RELERR,ASSERR, IFLAG,RES, 
JIJORK,  I WORK) 

IF  ( IFLAG.LE.3)  GO  TO  ISO 
IF  (IFLAG.EQ.9)  GO  TO  999 

FOUND  THE  MAX  CURRENT 


300  CONTINUE 

HVn(V-VLA5T) 

NHAX= IF IX(ALOG ( (V-VLAST) /(0 . 0 1 ) ) /ALOG (2 . 0) ) + 1 
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DO  A BINARY  SEARCH  FOR  THE  LIMITING  CURRENT 

II-  V'A  r.  "TINS  A D lU-NOT-CONVERGE  RETURN  CODE  V IS 
ASSUMED  TO  BE  LARGER  THAN  THE  LIMITING  CURRENT 

DO  500  N = 1 , NMAX 
DO  350  1=1, NEND 
D I5HAXC I ) =GAM'J- 1 . f) 

350  CONTINUE 

IF  C 1FL0G.LE.3)  GO  TO  =100 
DO  380  1=1, MEND 
C( I ) =GAM0 
380  CONTINUE 
400  CONTINUE 

FOR  THE  BINARY  SEARCH  ROUTINE  A REDUCTION  OF  10  ORDERS  OF 
MAGNITUDE  IS  CLOSE  ENOUGH 

HV=HV/2 . 0 

IF  ( IFLAG.LE.3)  VLAST-V 
IF  (IFLAG.LE.3)  V=V+HV 
IF  (IFLHG.GT.3)  V=V-HV 
I FLAG =1 

CALL  ON (F0FRZ,NR*NZ,C,M6AND,DISMAX, RELEER, ABSERR, I FLAG, RES, 
11J0RK,  I DORK) 

500  CONTINUE 

INSURE  THAT  THE  LAST  GOOD  CURRENT  IS  RETURNED  TO  THE  CALLER 
AND  RETURN 

600  CONTINUE 

IF  ( IFLAG.LE.2)  GO  TO  999 
IF  ( I FLAG . EC!  .3)  GO  TC  700 
V=VLAST 
NEND=NR*NZ 
DO  650  1=1, NEND 
DISMAX(I)=GAI-k3-I.0 
C ( I ) =GAM3 
650  CONTINUE 
IFL.AG=  1 
700  CONTINUE 

CALL  ON CFCFRZ, NR*NZ, C.MSAND, DISMAX, RELERR, A3SERR, 1FLAG, RES. 
ll,rORK,  I WORK) 

IF  ( I FLAG.  El).  3)  GO  TO  700 
399  CONTINUE 
RETURN 
END 
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SUBROUTINE  F0FRZ(C,F) 

comoN/Dflin/ra,  rh  , zn  , nr,  nrl  i , nrl2  , nz  , v,  cams 
COmON/RflDIUS/RdOl) 

C0m0N/LEN5TU/Z(  101) 

C0i-li-!0N/SUB2/ JFLAG,  H,  H,  K,  G C 7,  ?) 

REAL  K 

DIl-lENSION  C(  1C00) ,F(  1000) 

11=7 

DO  200  1=1, NR 
DO  100  J=1,HZ 

CHECK  TO  INSURE  THAT  EACH  VALUE  OF  GAMMA  IS  GREATER  THAN  1.0 
IF  (C(J+(I-1)*NZ).LE. 1.0)  GO  TO  900 


INETGRATE  OVER  QUADRANT  2 

F ( J-K  I- 1 ) *NZ) =RZ INT2 ( I, J, C) 

IF  (IFLAG.NE.U)  GO  TO  900 

INTEGRATE  OVER  QUADRANT  1 UNLESS  Z.GT.L/2 

IF  (J.LT.NZ)  F(J+( 1-1 )*NZ) =F(  J-K i-l)*NZ)+RZINTl ( I, J, C) 
IF  (IFLAG.NE.O)  GO  TO  900 

INTEGRATE  OVER  QUADRANT  3 UNLESS  R IS  TOO  SHALL 

IF  (I.GT.l)  F(J+(I-1)*NZ)=FCJ+(I-1)*NZ)+RZINT3(I,J,C) 
IF  (IFLAG.NE.O)  GO  TO  900 

INTEGRATE  OVER  QUADRANT  4 UNLESS  


IF  CU.GT.  1)  .AND.  (J.LT.NZ)) 

1 F ( J K I - J ) *NZ ) ~F  (J  f ( I~1 ) *NZ ) -f  RZ I NT4  (I , J , C ) 

IF  ( IFLAG.NE.U)  GO  TO  900 
100  CONTINUE 
200  CONTINUE 
RETURN 

NEGATIVE  SQUARE  ROOT////ASSUHE  V IS  GREATER  THAN  V LIMIT 

900  CONTINUE 
IEND=NR*NZ 
DO  910  1 = 1,  IEND 
F ( I ) = 1 . 0 
910  CONTINUE 
RETURN 
END 
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REAL  FUNCTION  R2INT1 ( I, J,C) 

CUl-TON/DATA/RO , RN . ZN , HR.  HRL 1 , NRL2 , NZ,  V,  GAM3 

common/radius/rubi) 

COMMON/LENGTH/Z  (10 1 ) 

COmON/SUNT/IF  LRG.fl.  H»K,G(?\7) 

DIMENSION  C(100O) 

REAL  K 
C 

C QUADRANT  1 INTEGRATOR 

C 

RZ INTI =0.0 

H=(R( I+1)-R(I) ) /FLOAT (IT- 1 ) 

K=(Z(J+2)  -2  C J+l ) ) /FLOAT  (11- 1 ) 

DR1=-1.0AR(J  + 1)-R(I)) 

DKI I =-DRI 

DZJ=- l . 0/(Z( J+2)-Z( J+l ) ) 

D2J)=-DZJ 

DO  200  L = 1 , 1*1 

RA=R ( I ) +FLOAT (L- 1 ) >:*H 

Rl=(RCI+l)-RA)/(R(Ivl)-R(I)) 

RI 1 E (Rfi-RC I))/(R(I+1)-R(I)) 

DO  100  N = 1 , 1*1 

ZArZ ( J+l ) +FLUAT  (N- 1 ) *l< 

ZJ=(Z(J+2)-ZA)/(Z(J+2)-Z(J+D) 

ZJl=(ZA-Z(J+l))/(Z(J+2)-Z(J+l)) 

Cl'C(J+(I-l)*tlZ) 

C2=C(J+HCI-1)*NZ) 

c'5=Gnm 

IF  (I.LT.NR)  C3=C(J+I*NZ) 

04=001*10 

IF  (I.LT.NR)  C4=C(J+1+I*NZ) 

F 1 = (ZJ**2)  *(PR  1**2)  +(R  1**2)  *(DZJ**2) 

F2=  (ZJX'ZJ  1 ) *<DR  I**2)+(RI**2)  *(DZJ*PZJ  1 ) 

F3=  ( ZJ“v"!*2 ) * ( DR  I *PR  1 1 ) + ( K 1 *E 1 1 ) * ( DZ  J **2 ) 

F4= ( Z J 1 *2 J ) *< DR I *PR II ) + ( R 11 *R 1 ) * ( PZ J 1 *DZJ ) 

G ( L , N)  = ( C H F 1 +C2+F2-I  C3*F3+C4*F4 ) *RA 
IF  (Cl. LT.MRL1) . OR. ( I .GE.NRL2) ) GO  TO  103 
GAM=Cl+RI*ZJ-tC2*RI*ZJl+C3*RIl*ZJ+C4*RI  1*2J1 
IF  (GAM. LE. 1.0)  GO  TO  AGO 

RAD  IC=  (4.  U*V*RA*GAiD  A (S0RT(GAM**2- 1.0))  *(RN**2-RG#*2) ) 
G(L,N)=GCL,N)+RADIC*RI*ZJ 
100  CONTINUE 
200  CONTINUE 

C 

C DO  A TRAPAZOIDAL  RULE  INTEGRAL 

C 

CALL  INTGLCRZINT1) 

I FLAG *0 
RETURN 

C 

C NEGATIVE  OR  ZERO  SQUARE  ROOT/AWNOTIFY  CALLING  PROGRAM 

C 

9G0  CONTINUE 
IFLAG=1 
RETURN 
END 
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REAL  FUNCTION  R2INT2(I.J,C) 

COMMON/DATA/RO, RN, ZN.NR. NRL  l , HRL2.NZ. V,  GAM3 
COMMON/'LENGTH/'ZC  101) 

COMMOIMiADIUS/RClOl) 

C0MM0N/’SUB2/iFLAC,M,H.K,G(?\7) 

REAL  K 

QUADRANT  2 INTEGRATOR 

DIMENSION  CC10O0) 

RZINT2=0.0 

H=(R(I-H)-R(I))  /FLOAT  (II- 1 ) 

K*>(Z(J-H)-Z(J))/FLOAT(M-l) 

DRI=-1.0AR(I  + 1)-R(I)) 

DR  1 1=-DR1 

DZJ-l.OAZU+D-ZCJ)) 

DZJ 1 =-DZJ 
DO  200  L = 1 , M 
RA-R(  I)+FL0AT(L-1 )*H 
RI=(R(I+l)-RA)/(R(i  H)-R(I)) 

R 1 1 = CRA-R ( I) )/(R( I + l)-R( I) ) 

DO  100  N = 1 , M 

zn=z(j)  tfloatcn-dn'K 

ZJ=(ZA-Z<J))/(Z(J+1)-ZCJ)) 

zji  = czcj+n-zn)/(zcj+i)-zij)) 

Cl -GAMS 

IF  (J.GT.l)  C1=C(J-  1+(I-1)*NZ) 

C2=C(J+CI-1)*NZ) 

C3-GAM3 

IF  ((J.GT.l). AND.C1.LT. NR))  C3-C(J-1+I*NZ) 

C4=GAM'J 

IF  (I.LT.NR)  C4-C(J+I*NZ) 

F1"(ZJ1  ZZ  J ) * (DR ! **2 ) + ( R I **2)  * <PZJ1*9ZJ) 
F2-(ZJ>R*2)*(DUI#*2)+(Rl:.-:-::2)  *(DZJ*#2) 

F3  - ( Z J 3 *Z  J ) * ( DR  I 1 *DR  1 ) + ( R II  *R  I)  * ( DZJ  1 *DZJ ) 
F4-CZJ**2)*(DRI1  T)RI)+(RI  1*RI)*(PZJ;:*2) 

G(L,N)  --(Cl *F  1 -i-C2*F2+Co*F3H C4*F 4 ) *RA 
IF  ( ( ! .LT.MRL 1 ) . OR. ( I . GE.NRL2) ) GOTO  100 
GAI1=Cl>CRI#2Jl+C2>l;RI-:€J+C3s:RI  J::<ZJ1+C4*RI  1*ZJ 
IF  (GAM. LE.  1.0)  GO  TO  .900 

RAD  1 C = (4. 0*V*RA*GAM) A (SORT (GP.IN-i-.2- 1.0)) *(RN**2~R0**2) ) 
G(L,N)=G(L,N)+RAPIC*R1*ZJ 
100  CONTINUE 
200  CONTINUE 

DO  A TRAPAZOIDAL  RULE  INTEGRAL 

CALL  INTGL(RZINT2) 

I FLAG  *=0 
RETURN 

NEGATIVE  OR  ZERO  SQUARE  ROOT/////NOTIFV  CALLING  PROGRAM 

900  CONTINUE 
IFLAGal 
RETURN 
END 
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REAL  FUNCTION  RZINT3(I, J,C) 

COMNON/DOTfl/Ra,  RN,  ZN, NR,  NRL 1 , NRL2,  NZ,  V,  GAMS 
COmON/LENGTH/ZUOl) 

COl'iMON-'ROD  I US/R  (101) 

C0m0N/SUB2/'IFLRG.ILH,K.G(?,7) 

REAL  K 

GUADRANT  3 INTEGRATOR 

DIMENSION  CC1GC0) 

RZINT3=0.Q 

H=(R(I)-R(I-i))/FL0AT(I1-l) 

K=(Z(J+l)-Z(J))/FLOflT(N-n 
DRI  = 1.0/(R(I)-R(I-D) 

DR  1 1=-DRI 

DZJ"1.0/'(ZCJ+1)-Z(J)> 

DZJ 1 =-DZJ 

do  2ue  l=i,h 

RA=R ( I - 1) +FLOAT (L- 1 ) *H 
R I - (RA-R ( I - 1 ) ) /(R ( I ) -R ( I- 1 ) ) 

R 1 1 = (R ( I) -RA)/(R( I)-R(I-l)) 

DO  IOC  N=  1 , N 

ZA=Z(J)+FL0ATCN-1)*K 

Z J = ( ZA-Z  ( J ) ) / ( Z ( J+ 1 ) -Z  ( J ) ) 

zji  = (Z(j+n-ZA)/cz(j+i)-z(j)) 

Cl-GAHO 

IF  CJ.GT.l)  Cl“C(J-l+( 1-2) *NZ) 

C2«C(J-KI-2)*NZ) 

C3=GAH3 

IF  CJ.GT.I)  C3=C(J-1+(I-1)*NZ) 

C4=C(J+(I-n*NZ) 

F 1 = (ZJK'-ZJ ) *(PR  1 1-DR  I ) -KR 1 1 *R I ) >:  (DZJ  1 *DZJ) 

F2=(7J**2) *CDR  I ! X-DR I ) + CR I 1*R  I ) *CDZJ**2) 

F3*  CZJ 1 *ZJ  l * ( DR  I »*2 ) -HR  I **2 ) * ( DZJ 1 *DZJ ) 

F4«(ZJsM2)*(DR 1**2)  •)  (RI**2)  *( DZJx  *2) 

G(L,N)  “(Clx-F  l+C2;!F2+C3'i<F3+r.4+F4)  *RA 
IF  ((I.LE.NRL1)  .OR.  ( I .GT.IIRL2) ) GO  TO  100 
GAN=C  I X R 1 1 *ZJ  1 +C2>;  R 1 1 *Z  J+C3*R  I •"Z  J 1 +C4*R 1 -NZ  J 
IF  (GAM.LE. 1.0)  GO  TO  903 

RADIC- (4. 0*V*RA*GAM)/( (SCRT(GnM**2- 1.0)) *(RN**2  R0**2) ) 
G(L,N)=G(L,N)+RADIC*Rl*ZJ 
100  CONTINUE 
200  CONTINUE 

DO  A TKAPAZOIDAL  RULE  INTEGRAL 

CALL  INTGLCRZINT3) 

IFLAGe0 

RETURN 

NEGATIVE  OR  ZERO  SQUARE  K’OOT/AV/NOTIFY  CALLING  FROGRAM 

SCO  CONTINUE 
I FLAG =1 
RETURN 
END 
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REAL  FUNCTION  RZINT4C I, J,C) 

CONraN/PATfVRO,  RN,  ZN,NR,NRL  1 , NRL2,  NZ,  V,  GAT13 
COI'TON/RAD  I IJS/R  (101) 

COIlI-lON/LCNGTH/ZUOl) 

C0!TON/SUB2/'IrLAG,N,H.K,G(7.7) 

REAL.  K 

QUADRANT  4 INTEGRATOR 


DIMENSION  C( 1030) 

RZIN73=4.0 

H=  (R  ( 1 ) -R  ( I- 1 ) ) /FLOAT  (II- 1 ) 

K“(Z(J+2)-ZCJ+l))/FL0AT(M-l) 

DRI=1.0/(R(  I)-R(I-l)) 

DR  1 1 B-PR I 

PZJ  — 1.0/(ZCJ+2)-Z(J+l)) 

DZJl'-DZJ 

DO  200  L = 1 , M 

RA=R ( I- 1 ) +FLOAT (L-  1)*H 

R I = (RA-R ( I- 1 ) ) /(R ( I) -R ( I- I ) ) 

RI1 =(R( I )-RA)/(R( I)-R( 1-1) ) 

DO  100  M = Ml 

ZA“Z ( J+l ) TFLOAT (N- 1 ) *K 

ZJ=CZ(J+2)-ZP)/(Z(J+2)-Z(J+l)) 

ZJl“(ZA-Z(J-:-l))/(Z(J+2)-Z(J+l)) 

ClBC(J+(I~2)*NZ) 

C2=C( J+l+(  I -2)  'JNZ) 

C3bC(J+(  I-1)!!-NZ) 

C4---CU+1-H  I-1)*NZ) 

F 1 = (ZJ**2)  *(DR  11*021)  -HR  I MR  I ) *(DZJ**2) 

F2  = ( Z J 1 vZJ ) * ( Dim  >:•!*!’  n + (R 1 1 *R  I ) * C DZJ 1 *DZJ ) 

F3»  (ZJ#*2)  >:*<DR  I>‘~*2) -HR  I**2)  *<DZJ**2) 

F4-CZJ)  *ZJ ) * ( DR  I ;!- : 2)  H R I **2 ) #( DZ J 1 *PZJ) 

G (L,  N)  ••  (C 1 ‘ F 1 +C?:'  F 2 +C 3 ■ F 3 +C- 4*F 4 ) *RA 
IF  ( ( 1 .LF.  .URL  1)  . OR.  (I  .GT.NRL2) ) GO  TO  100 
GAi'l=C  1*R11  :-.-2J+C?*R  1 1 *2J  1 +C3F  R I *ZJ  +C4*R  I *ZJ  1 
IF  (GAI1.LE.  1.0)  GO  Til  000 

RAD  ICB  (4. 0*V:2A-TGAM)  /( CS0RT(GAI1**2- 1.0)  *(RN**2-R0**2) ) 
G(L,N)-G(L, N) TRAD  I C*R 1 *ZJ 
IPO  CONTINUE 
200  CONTINUE 

DO  A TRAFAZOIDAL  RULE  INTEGRAL 

CALL  INTEL (RZIHT4) 

IFLAG-0 

RETURN 

NEGATIVE  OR  ZERO  SQUARE  ROOT/////NUTIFY  CALLING  PROGRAM 

900  CONTINUE 
I FLAG «1 
RETURN 
END 
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SUBROUTINE  INTEL (RESULT) 

DO  P TRAPAZOIPAL  INTEGRAL  APPROXIMATION  ON  THE  DATA  IN  G 

C0MM0N/SUA2/IFLAG,  II,  H,  K,G  (.7,  7) 

REAL  K 

RESULT- U.O-M.QWGU,  1)+G(  1,KHG(N,  1)+G(M.M>) 

JENDCM-1 
DO  100  Jv2. JEND 

RESULT “RESULT  + ( 1 . 0/2 . 0) *<G ( 1 , J) +G ( J , 1 ) +G (M.  J) +G ( J,  N) ) 

J0O  CONTINUE 
1END=JEND 
DO  300  I =2, IENP 
DO  200  J-2,JENP 
RESULTrRESULT+G( I, J) 

200  CONTINUE 
300  CONTINUE 

RESULT-H«*RESULT 

RETURN 

END 
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